The Forward Model

To find the relation between a dosing regimen and an observed effect, the relation between the dosing regimen and time-course of the drug concentration inside the body must initially be determined. Pharmacokinetic (PK) models determine this relation through mathematical models of the transport of the drug around the body. Once this time course is found, It can be inputted into a Pharmacodynamic (PD) model, which translates the drug concentration into a measurable PD effect

The PK Model

The empirical PK model that is most commonly used separates the body into a number of compartments in which the drug can reside, see Figure \ref{fig:generalPK}. There is the central compartment, usually representing the blood, and a number of peripheral compartments, representing different groups of tissues. The drug in the central compartment can be transferred to the peripheral compartments and vice versa, and it is also cleared from the central compartment. Additionally, there may be a dosing compartment depending on the method of dose administration. The $n$-compartment PK model is given by the following system of ordinary differential equations (ODEs), \begin{align*} \dot{D} =& - Q_d\left(\frac{D}{V_d}\right), \\ \dot{C} =& Q_d\left(\frac{D}{V_d}\right) - CL\left(\frac{C}{V_c}\right) - \sum_{i=1}^{n-1}Q_{i}\left(\frac{P_i}{V_i}, \frac{C}{V_c}\right),\\ \dot{P_i} =& Q_{i}\left(\frac{P_i}{V_i}, \frac{C}{V_c}\right), && 1\le i \le n, \end{align*} where $D$, $C$ and $P_i$ are the masses of drug in the dosing, central and $i$-th peripheral compartments, respectively, $V_d$, $V_c$ and $V_i$ are the volumes of the dosing, central and $i$-th peripheral compartments, respectively, $Q_{i}$ is a function describing the flux of the drug into the $i$-th peripheral compartment, $Q_d$ is a non-negative function describing the flux of the drug into the central compartment from the dosing compartment, and $CL$ is a non-negative function describing the loss of the drug from the central compartment through clearance. The peripheral compartments have initial conditions $P_i(0)=0$, the dosing compartment has initial condition $D(0)=d$, the dose amount, or $D(0) = 0$, depending on the method of dose administration. Similarly, the central compartment has initial condition $C(0)=0$ or $C(0) = d$ depending on the method of dose administration \cite{mould2012, Mould2013, nielsen2013}.

Figure 1: Diagram of the general $n$-compartment PK model. The drug enters the body either through a dosing compartment, or straight into the blood stream, represented as the central compartment. It then moves in and out of different tissues, represented as peripheral compartments, before being cleared from the central compartment.

This gives a family of models which has various resulting concentration-time curves. To determine the model that should be used for simulation there are three decisions to make:

The dosing compartment is an easy decision to make as it depends on the form of drug administration that is being modelled. For example, if the drug is administered via IV, then no dosing compartment is required, but if the drug is administered subcutaneously then the dosing compartment is required. The number of compartments and the form of the transfer rate are decided based on the shape of the output that is needed. This can be best determined through the time vs. drug concentration curve plotted on a log scale (see figure 2). If one compartment is used, this plot is a straight line (due to there being a single exponential term in the solution). If two compartments are used then two straight sections can be seen in the plot (due to two exponential terms in the solution). And for $n$-compartments, there can be up to $n$ straight sections in the plot. A model of $n$ compartments can produce a curve with less than $n$ straight sections however these dynamics can also be reproduced with a simpler model. For deciding the form of the rates used, the effect of different doses need to be plotted (see figure 3). If all rates are linear, i.e.

\begin{align*} Q_d\left(\frac{D}{V_d}\right) &= k_d\frac{D}{V_d}, \\ CL\left(\frac{C}{V_c}\right) &= k_0\frac{C}{V_c},\\ Q_{i}\left(\frac{P_i}{V_i}, \frac{C}{V_c}\right) & = k_i(\frac{C}{V_c} - \frac{P_i}{V_i}), \end{align*}

then the traces of the different doses would be parallel. However, if [...].

Figure 2: Simulations of one compartment (left), two compartment (middle) and three compartment (right) PK models with linear rates of transfer between compartments and no dosing compartment. Y-axis is on a log-scale

In practice, the model choice comes from the data that is collected. For example, the data shown in figure 4, is from IV bolus of docetaxel in rats, clearly shows two straight sections and the dose-normalised curves mostly line up. Thus the most appropriate model for this drug would be a linear two-compartment model with no dosing compartment.

Figure 3: Experimental data of rats given different levels of docetaxel. Observations (Docetaxel concentration, ng/mL) have been normalised by dose amount and then averaged over the rats in each dose group.

The PD Model

To translates this concentration of drug in the central compartment, $C_\mathrm{conc}(t)$, into a direct drug effect, $E_\textrm{drug}$, commonly the PD relation is either a linear relation, $E_\textrm{drug} = s C_\mathrm{conc}(t)$, where $s>0$ is the slope of the linear effect; or a non-linear $E_\textrm{max}$ relation, \begin{equation*} E_\mathrm{drug}(t) = \frac{E_\textrm{max}C_\mathrm{drug}(t)}{C_\mathrm{drug}(t) + EC_{50}}. \end{equation*} where $E_\textrm{max} > 0$ is the maximum of the effect and $EC_{50}$ is the concentration of drug at which the effect is at 50\% of $E_\textrm{max}$.

Unfortunately the direct effect is not always observable and an indirect response model, which relates drug effect to drug response, is required. This is the case for myleotoxicity where the drug effect, the level of inhibition of proliferation, is not directly measurable but the response, the number of blood cells in circulation, is measurable. The model presented by Dayneka \textit{et al.} \cite{dayneka1993comparison} captures this type of indirect response. In this paper, the authors present four different models of action a drug could take on the response variable, $R$, the variable of interest to the researcher. In the myelotoxicity case, $R$ is the concentration of the blood cell that is being modelled. These models are of a general biological system where the response variable is produced at a constant rate of $k_{\textrm{in}} > 0$ and eliminated at a rate $k_{\textrm{out}} > 0$. The drug could affect either the production or the loss of $R$; and it is either inhibitory or stimulatory. For myelotoxicity, the production of blood cells is inhibited by the drug, and so the differential equation for the response variable is $$\dot{R} = k_{\textrm{in}}(1-E_\textrm{drug}(t)) - k_{\textrm{out}}R ,$$ where $E_\textrm{drug}(t) < 1$ is the inhibitory drug effect. The initial condition, $R(0) = R_0$ is such that the system begins at the drug-free steady state, i.e. $k_{\textrm{in}} = k_{\textrm{out}} R_0$.

However, there are some features seen in myelotoxicity data that this general model cannot capture, such as a delay between drug administration and the observed response of a drop in blood cell concentrations, as can be seen in Figure \ref{fig:ExampleData}. This delay comes from the conversion of proliferating cells into circulating cells which requires time for the cells to mature and enter circulation. As the proliferating cells are affected by the drug but the circulating cells are observed, this delay is seen in the data. Measuring the response in bone marrow directly in patients is not possible, thus a more specific model is required to capture this delay.

One such model has been presented by Zamboni \textit{et al.} \cite{zamboni2001pharmacodynamic} who extended the model from one to three compartments (see Figure \ref{fig:allPD}B). One compartment represents the proliferating cells, $P$, which the drug affects the production of; another compartment represents the measurable cell count in circulation around the blood, $C_{\mathrm{cell}}$; inbetween these is a transit compartment to represent the time delay. This leads to three differential equations: \begin{align*} \dot{P} & = k_{\textrm{in}}E_\mathrm{drug}(t) - k_{\textrm{tr}}P, \\ \dot{T} & = k_{\textrm{tr}}P - k_{\textrm{tr}}T, \\ \dot{C_{\mathrm{cell}}} & = k_{\textrm{tr}}T - k_{\textrm{out}}C_{\mathrm{cell}}, \end{align*} where $k_{\textrm{tr}} > 0$ is the rate of transfer between the compartments. The initial conditions were not provided in the paper.

Noise Model